/*!
 ***************************************************************************
 * \file   heat2.cpp 
 * Solving the Heat difusion equation in 2D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Sun Apr 13 19:48:52 BST 2008 ]
 ***************************************************************************
 * This code solves the heat diffusion equation in two dimensions: 
 * \f$ \frac{\partial T}{\partial t} = \nabla^2 T\f$. 
 * \par
 * This equation is solved in a unit square 
 * \f$ x, y \in [0,1] \f$. The boundary conditions are as shown in
 * the next figure:
 *
 * \image html LaplaceSquare.png "Unit square in 2D" width=10cm
 * \image latex LaplaceSquare.eps "Unit square in 2D" width=10cm
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./heat2  \endverbatim
 * The results are stored in \c Data2 .
 ***************************************************************************
 */
#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"

using namespace Tuna;

// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

int main() 
{
    double length_x, length_y;
    int num_nodes_x, num_nodes_y;
    int num_vols_x, num_vols_y;
    double tolerance, dt;    
    int steps, tdma_iter, frequency;

// Reading data from the input file
    std::ifstream input_cin ("input2"); 
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y
	      >> dt
	      >> steps
	      >> tolerance
	      >> tdma_iter;
    input_cin.close();

    StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    mesh.print();
    ScalarField2D T ( mesh.getExtentVolumes() );
    ScalarField2D Tn( mesh.getExtentNodes()   );

    Range I ( Tn.lbound(firstDim), Tn.ubound(firstDim) - 1);
    firstIndex i;
    Tn(I, Tn.ubound(secondDim)) = 10 * sin (PI * dx * i);
    NumUtils::interpolateToVolumes(T, Tn);


#ifdef WITH_DX
//    InOut::writeToFile_DX(Tn, 0, "./Data2/temp.", dx, dy);
    InOut::writeToFile(Tn, 0, "./Data2/temp.", dx, dy);
#endif

#ifdef WITH_GNUPLOT
    InOut::writeToFile(Tn, 0, "./Data2/temp.", dx, dy);
    GNUplot plotter;
    plotter("set grid");
    plotter("splot \"./Data2/temp.00000\" title \"Initial Temperature\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif

    SparseMatrix< Diagonal <double, 2> > A(num_nodes_x, num_nodes_y); 

    ScalarField2D               b(num_nodes_x, num_nodes_y);
    
    ScalarEquation<CDS<double, 2> > heat(T, A, b, mesh.getDeltas());
    heat.setDeltaTime(dt);
    heat.setGamma(1.0);
    heat.setDirichlet(TOP_WALL);
    heat.setDirichlet(LEFT_WALL);
    heat.setDirichlet(RIGHT_WALL);
    heat.setDirichlet(BOTTOM_WALL);
    heat.print();

    ofstream error_fp("./Data2/error"), residual_fp("./Data2/residual");
    int iteration = 0;
    double error = 1.0, residual;
   
    while ( (error > tolerance) && (iteration < steps) ) {
	heat.calcCoefficients();
	Solver::TDMA2DY(heat, tolerance, tdma_iter);
	residual = heat.calcResidual();
	error = NumUtils::calcErrorL1(T, heat.phi);

	error_fp << iteration << "\t" << error << endl;
	residual_fp << iteration << "\t" << residual << endl;

     	cout << "\n Iterations " << ++iteration 
	     << " Residual = " << residual 
	     << " Error = " << error;

#ifdef WITH_DX
    NumUtils::interpolateToNodes(Tn, T);
//    InOut::writeToFile_DX(Tn, iteration, "./Data2/temp.", dx, dy);
    InOut::writeToFile(Tn, iteration, "./Data2/temp.", dx, dy);
#endif

	heat.update();
    }

#ifdef WITH_GNUPLOT
    NumUtils::interpolateToNodes(Tn, T);
    InOut::writeToFile(Tn, 1, "./Data2/temp.", dx, dy);
    std::cout << "\n >---> Press <1> and then <enter> to plot the solution \n";
    std::cin >> pause;	
    plotter("set contour");
    plotter("set cntrparam level 20");
    plotter("set view 0,0,1.5,1");
    plotter("set size 0.73,1");
    plotter("unset surface");
    plotter("splot \"Data2/temp.00001\" title \"Final solution\" w l");
    std::cin >> pause;
#endif

#ifdef WITH_DX
	std::cout << "\n\n NOTE: The data were written in the [./Data2] directory in DX format";
#endif

    return 0;
}




